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Abstract 


A summary of research conducted during the first year is presented. The research objectives were 
sought by conducting two tasks: (1) investigation of probabilistic design techniques for reliability- 
based design of composite sandwich panels, and (2) examination of strain energy density failure 
criterion in conjunction with response surface methodology for global-local design of damage 
tolerant helicopter fuselage structures. This report primarily discusses the efforts surrounding the 
first task and provides a discussion of some preliminary work involving the second task. 

1. Reliability-Based Design of Composite Sandwich Panels 

A reliability-based design of a rectangular sandwich plate with anisotropic face sheets and edges 
elastically restrained against rotation is investigated. The plate, as shown in Fig. 1, is optimized for 
minimum weight subject to a reliability constraint against axial and shear buckling. The influence of 
reliability requirement on the optimal design is evident by comparing the results of the probabilistic 
design with those of a deterministic design. 

The plate is modeled based on the general small-deflection theory of sandwich plates with the mean 
values of the in-plane buckling loads determined using the Rayleigh-Ritz method according to the 
procedure discussed by Marcellier and Rais-Rohani. The plate design is optimized using the 
modified method of feasible directions in DOT 2 software. 



Figure 1. General description of the sandwich plate model 
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1.1 Reliability Analysis 


The reliability is defined here as the probability that the panel will not buckle under the applied load. 
Mathematically, this statement can be expressed as 

R e = P(N cr >N a ) CD 

In this context, the performance function that describes plate reliability is given as 

Z = N cr -N a =g(X) ( 2 ) 

where X is the vector of basic random variables affecting plate reliability, and g(.) is a function that 
describes the relationship between the performance function and the basic random variables. Plate 
failure is defined by Z < 0; plate survival is defined by Z > 0; and the limit state is defined by Z = 0. 

In advanced second-moment method (ASM) for calculating plate reliability, the basic random 
variables are transformed into a normalized coordinate system according to the relation 

Xj = (x, -/r-J/ a- for /=1,2, ..., n (3) 

The resulting normalized random variables have mean of zero and standard deviation of one (i.e., 
= Oand Oy = 1). Based on ASM the shortest distance, in the reduced coordinate system, 

measured from the origin to the limit state (Z = 0) or failure surface defines the reliability index (3 
also known as the Hasofer-and-Lind (HL) index. 

The iterative procedure described by Ayyub and McCuen is used here to calculate /i. The plate 
reliability is then calculated using the relationship between ft and the probability of failure given by 

p f = \-&(P) = l-R e (4) 

The decision on what reliability method to use followed an initial comparison between ASM and the 
first order reliability method (FORM). Although the latter is easier to implement, it was found to 
underestimate plate reliability mainly due to the non-linearity of the performance function. Hence, 
only ASM was used in reliability-based optimization. 

1.2 Effectivitv Analysis Based on Taguchi Method 

In order to identify the statistically significant parameters affecting the buckling response of a 
sandwich plate, a design of experiments (DOE) was set up based on the principles of Taguchi 
method. Eleven main factors at two levels and four two-factor interactions were considered in the 
DOE. Using the L 16 orthogonal array a total of sixteen experiments were conducted. The order of 
factors influencing the panel buckling strength is found to be h c , t, t*h c , K, t*K, G 12 , G zx , 0, G yz , v, 
t*0, 12,, E p (a/b)*K, and a/b. The description of each parameter and its mean value are given in 
Table 1. Among all factors the face sheet ply thickness and the core thickness are found to have the 
most significant influence on the buckling load. Details of this investigation including the robust 
design configuration for buckling load maximization can be found in ref. 4. 


2 



Table 1. Description of basic random variables 


Basic Random Variable 

Mean 

CoV, % 

Distribution 

axial load, N x 

2000 N/mm 

2.5 

Normal 

axial load, N w 

2000 N/mm 

2.5 

Normal 

plate length, a 

254 mm 

5.0 

Normal 

plate width, b 

254 mm 

5.0 

Normal 

edge rotational stiffness along 

10 

5.0 

Normal 

y = 0, b edges, K = bKJD 22 
Face Sheets: graohite-epoxy 

ply thickness 3 

t()> Us’ t-45* bo 

5.0 

Normal 

ply angle b 

0, 45, -45, 90 

5.0 

Normal 

Young’s modulus, E, 

229 GPa 

2.5 

Normal 

Young’s modulus, E 2 

13.35 GPa 

2.5 

Normal 

shear modulus, G n 

5.25 GPa 

2.5 

Normal 

Poisson’s ratio, v n 

0.315 

2.5 

Normal 

Core: aluminum honeycomb 

thickness 3 

K 

5.0 

Normal 

shear rigidity, G x , 

0.146 GPa 

2.5 

Normal 

shear rigidity, G„ 

0.0904 GPa 

2.5 

Normal 


a design variables 

b For the 0° ply a standard deviation of ±1° is used 


1.3 Probabilistic Design Optimization 

For the composite sandwich plate the probabilistic optimization problem is formulated as 
Min. ftX,Y) 

St. p (5) 

Y l <Y <Y U 

where the objective function /is the plate weight, F is the vector of design variables, representing a 
subset of the random variable vector X, defined in Table 1, and P min is the required minimum 
reliability index which in this case is set to 3.09 for a reliability of 0.999. All random variables are 
assumed to have normal distributions with specified means and assumed coefficients of variation 
(CoV) given in Table 1 . 

Each face sheet is an unsymmetric quasi-isotropic laminate with four plies made of graphite-epoxy 
material. The sandwich plate is symmetric about its midplane surface. The mean values for ply and 
core thicknesses are treated as the only design variables in the optimization process. In the 
probabilistic design problem the performance function in Eq. (2) is approximated as 

H z = H(N cr )-^(N a ) ( 6 ) 

with mean value of the buckling load obtained from the deterministic procedure described in ref. 1 . 
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1.4 Deterministic Design Optimization 

For the composite sandwich plate the deterministic optimization problem is formulated as 
Min. J{X,Y) 

S.t. N cr > FS (N a ) ( 7 ) 

Y 1 <Y < Y u 

where the objective function /is the plate weight, and FS is the factor of safety against buckling. 
The design variables in this case are assumed to be deterministic. The design reliability is 
calculated for the optimal design as a post-optimization step. 

1.5 Results and Discussion 

The deterministic and probabilistic optimization problems are solved each for two different loading 
conditions (N x and A/). In each case the bounds on design variables are: 0.0254 mm < t ply < 7.62 
mm and 0.254 mm <\ < 16.5 1 mm. The face sheet material has a density of 1600 kg/m 3 while the 
honeycomb core has an effective density of 27.1 kg/m 3 . 

To demonstrate the effect of rotational edge stiffness, three different plate boundary conditions are 
considered. In the first case the plate is assumed to be simply supported, i.e., K, = K 3 = K 2 = K 4 = 
0. In the second case the plate is simply supported at x = 0 and a edges but elastically restrained 
against rotation along the other two edges such that K, = K 3 — 0 and K 2 = K 4 = 10. In the last case, 
the rotational stiffness along y = 0, b edges is increased to infinity making those two edges 
clamped. For each case, the optimal design is obtained for three different aspect ratios {a/b= 0.5, 1.0, 
and 2.0). The initial design for all cases is chosen as t 0 = t 45 = t_ 45 = t 90 = 1.27 mm and h c =15.24 
mm. 

The results of optimization are given in Tables 2 and 3. In all cases the core thickness is pushed 
near its upper bound. Since thickness has a large influence on plate bending stiffness, it is 
reasonable to see the core thickness increase. This increase in thickness causes the face sheets to 
move farther apart, therefore, increasing the plate bending stiffness D Vj . 

However, a similar increase in thickness is not observed in the face sheet plies. This is because the 
density of face sheet plies is significantly higher than that of the core. Hence, increasing the face 
sheet thickness has a greater impact on the overall weight of the plate that is being minimized. 

The thickness is seen to increase only in plies that offer greater support against buckling. For 
example, the contribution of the 90° plies against axial buckling is far less than the others. 
Therefore, we see this ply thickness remain near its lower bound value in almost all optimal designs. 

The deterministic and probabilistic optimal designs for axial loading condition indicate the 0° plies 
to be the most dominant followed by 45° and 45° plies. For the positive shear loading case, 
however, the 45° plies are more dominant than the 0° plies. It is expected that for the negative shear 
loading the 45° plies would be the most dominant. There is also a greater balance between the 
±45° plies in the axial than in the shear buckling case for all aspect ratios and support conditions 
investigated. 
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Table 2. Optimal design of (07±45Y90YCore) s sandwich plates under uniaxial compression, N x 
(K, =K, = 0; AY = K ± = bK/D „) 


bk/D 22 

t () , mm 

t 45 , mm 

V 

t, 4S , mm 

1 mm 

h c , mm 

4 22/ 

Deterministic Design 

Probabilistic Design 


Reliability Weight, gm 

Reliability 

Weight, gm 

a/b = 0.5, Area = 6.4516 x 10 2 m 2 

0 

0.2750 

0.0527 

0.0537 

0.0254 

16.51 

0.4920 

112.84 




0.5096 

0.1683 

0.0254 

0.0254 

16.50 



0.9990 

179.24 

10 

0.2798 

0.0363 

0.0365 

0.0254 

16.51 

0.4920 

106.90 




0.4542 

0.0876 

0.0888 

0.0254 

16.51 



0.9990 

164.26 

oo 

0.2418 

0.0342 

0.0340 

0.0254 

16.51 

0.4920 

98.10 




0.3808 

0.0752 

0.0766 

0.0254 

16.51 



0.9990 

144.04 




a/b = 1 .0, Area = 6. 

4516 x 10 2 m 2 




0 

0.2594 

0.1584 

0.1613 

0.0254 

15.29 

0.4920 

151.49 




0.4237 

0.1562 

0.1635 

0.0271 

16.33 



0.9987 

187.60 

10 

0.3297 

0.0883 

0.0929 

0.0254 

15.30 

0.4880 

137.46 




0.4493 

0.0887 

0.0912 

0.0254 

16.48 



0.9987 

163.95 

oo 

0.2331 

0.0421 

0.0490 

0.0254 

16.36 

0.4980 

100.76 




0.3810 

0.0728 

0.0785 

0.0254 

16.51 



0.9990 

143.97 




a/b = 2.0, Area = 6. 

4516 x 10 2 m 2 




0 

0.2382 

0.1224 

0.1262 

0.0254 

16.30 

0.4880 

134.255 




0.4205 

0.2199 

0.2156 

0.0271 

16.08 



0.9991 

210.39 

10 

0.2578 

0.0819 

0.0869 

0.0254 

16.51 

0.4900 

121.09 




0.3750 

0.1656 

0.1759 

0.0254 

16.37 



0.9989 

181.75 

oo 

0.2148 

0.0694 

0.0678 

0.0254 

16.50 

0.4900 

106.76 




0.3404 

0.1129 

0.1169 

0.0254 

16.50 



0.9990 

151.77 


The deterministic results indicate the axial loading to be a more severe condition than shear loading 
for all boundary conditions and aspect ratios. A similar trend is also observed in the probabilistic 
results with a few exceptions. The weight of the simply-supported plate under shear loading is 
found to be higher than that under axial loading for alb = 1 and 2. 

For both loading conditions, with a few exceptions, increasing plate aspect ratio (for a fixed plate 
area) causes the weight to increase. On the other hand increasing the rotational rigidity along the y = 
0, b edges causes the weight to decrease due to an increase in plate’s resistance against buckling. 

With FS = 1.0 in Eq. (7), the reliability of each deterministic optimal design is found to be roughly 
50%. The tolerance specified for buckling constraint in the optimization analysis is set at 0.001, 
and that is why in most cases the reliability is slightly less than 0.5. The reliability of each 
probabilistic optimal design, by contrast, is found to be at or near 0.999. 

For comparison, the plate having bk/D 22 = 10 and a/b = 1 was optimized for uniaxial compression 
with FS = 1.25. The optimal weight is found to be 385.71 gm which is approximately 180% higher 
than the corresponding case in Table 2. It is evident that the 25% increase in the factor of safety 
results in a substantial increase in the weight. 

Additional details of this investigation can be found in ref. 5. Future investigation will focus on 
maximization of reliability index and factor of safety (for deterministic design). This approach will 
determine the conditions under which the reliability can be increased to its maximum for a specified 
plate weight. 
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Table 3. Optimal design of (07±457907Core), sandwich plates under uniform shear load, N iy 
(K, = K, = 0-K 1 = K a =b&D„) 

bk/D 22 t„, mm t 45 , mm t. 4S , mm t w „ mm h„ mm Deterministic Design Probabilistic Design 


Reliability Weight, pm Reliability Weight. 


a/b = 0.5, Area = 6. 4516 x 10 2 m 2 


0 

10 

oo 

0.0254 

0.0548 

0.0254 

0.0372 

0.0254 

0.0436 

0.0386 

0.0819 

0.0257 

0.0729 

0.0254 

0.0566 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

0.0254 

16.51 

16.51 

16.51 

16.14 

16.51 

16.51 

0.5000 

0.4920 

0.5199 

52.58 

49.92 

49.84 

0.9987 

0.9989 

0.9990 

67.58 

61.42 

60.04 

a/b = 1.0, Area = 6. 4516 x 10 2 m 2 

0 

0.0621 

0.2578 

0.0254 

0.0254 

16.51 

0.4841 

105.373 




0.3518 

1.0380 

0.0254 

0.0254 

16.51 



0.9990 

326.22 

10 

0.0403 

0.2066 

0.0254 

0.0254 

16.51 

0.4860 

90.320 




0.1073 

0.4303 

0.0254 

0.0254 

16.50 



0.9988 

150.29 

oo 

0.0358 

0.1855 

0.0254 

0.0254 

16.50 

0.4860 

85.000 




0.1079 

0.2945 

0.0254 

0.0761 

16.50 



0.9990 

132.86 

a/b = 2.0, Area = 6. 45 16 x 10 2 m 2 

0 

0.0254 

0.3817 

0.0254 

0.0254 

16.51 

0.4880 

123.39 




0.1545 

1.0867 

0.0254 

0.0657 

16.51 



0.9990 

303.88 

10 

0.0254 

0.3555 

0.0254 

0.0254 

16.51 

0.4920 

117.98 




0.1189 

0.7096 

0.0254 

0.0968 

16.51 



0.9989 

225.09 

oo 

0.0254 

0.2723 

0.0254 

0.0264 

16.51 

0.5000 

100.99 




0.0548 

0.5383 

0.0254 

0.0690 

16.51 



0.9990 

170.78 
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2. Global-Local Design of Damage Tolerant Structures 


The second task in this research investigation involves the development of an efficient methodology 
for global-local analysis and optimization of damage-tolerant fuselage structures. The main focus 
areas are: (1) the application of strain energy density failure criterion, and (2) the development of 
response surface models for calculation of panel buckling load for use in global optimization. 


The finite element method is used for the global analysis of the structure subject to a 
multidimensional loading condition as shown by the model in Fig. 2. The structure in Fig. 2 (a) 
represents the global model and the panel in Fig. 2 (b) represents the local model. The global 
model is assumed to be simply-supported at the ends with displacements restrained in radial and 
angular directions and free in the z direction. The local model is assumed to be simply-supported 
along all edges. The global structural model is to be optimized subject to constraints imposed on 
maximum strain energy density in each shell element, maximum stress in each bar element, and 
buckling of each panel. 



Figure 2. (a) Cylindrical fuselage section; (b) Stiffened panel 

2.1 Strain Energy Density Failure Criterion 

The strain energy density (SED) failure criterion, first proposed by Sih 6 , is different from other 
classical fracture mechanics methods in that it does not require the presence of an initial flaw of 
known size and location. This criterion is based on the fact that strain energy stored in a material 
can be divided into dilatational and distortional parts, with failure initiation being associated with the 
dilatational part. 

SED criterion identifies critical regions in the structure where the dilatational part of strain energy is 
the more dominant part. These regions are found to be susceptible to damage initiation and 
subsequent failure of the structure. 

SED criterion of failure can be stated by the following two hypotheses 7 : 

1) Failure is assumed to coincide at locations where the local minima of strain energy density 
function (dW/dV) mjn is at its maximum. 

2) Failure initiation occurs when (dW/dV) mjn reaches its critical value governed by 


7 



dW 

~dV 


= const. 


(8) 




r \ 


r 2 


r 

( 


Where S is the critical strain energy density factor that serves as a measure of the fracture 
toughness value of the material and is found using 


S, = 


(l-2v)(l + v) 

2tt(1-v 2 ) ,c 


(9) 


where G Ic is the critical energy release rate of the material and v is the Poisson’s ratio. The 
unstable fracture begins when the ligament size reaches the critical value designated by r c . G lc is 
related to the mode I stress intensity factor according to 8 


G /( =d-v 2 )^ 


( 10 ) 


where K ]c is the critical stress intensity factor for mode I fracture of the material. 

The SED criterion has been used for failure prediction in both isotropic and composite materials. 
In the case of laminated composite structures, it is used to obtain information on the onset of 
delamination 7 . In that case r represents the distance from the delamination front. The onset of 
unstable delamination is assumed to coincide with (dW/dV) mjn reaching its critical value. 


Prior to the application of this criterion to the composite fuselage structure, it was tested on a 
20" x 10 " x 0.5" aluminum rectangular plate with a 0.1" circular hole at the center. The plate is 
considered to be simply supported along exterior edges, and is under the action of unifonn axial 
tensile force. Due to biaxial symmetry, one quarter of the plate, as shown in Fig. 3, is modeled 
for computational analysis. The model contains 120 four-noded quadrilateral shell elements with 
882 degrees of freedom. 



> 

x 


Figure 3. Finite element model of the plate with central circular hole 
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Symmetric boundary conditions are imposed along the axes of plate symmetry. The loaded edge 
is constrained against translation in the y and z directions while the unloaded edge is constrained 
only against translation in the z direction. 

To determine the location of failure initiation, the strain energy density contours are first 
generated. Figure 4 shows the SED contours near the hole corresponding to the axial compressive 
load of 250 lb. 





Index 


b 

n 

- □ 

- — 9 

... 8 

7 

6 


8. 237E-D4 
7. 553E-D4 
6.819E-04 
6.D36E-04 
5.352E-04 
4.618E-04 
3.884E-D4 
3. 150E-0 A 
2.416E-04 
1.682E-04 
9. 433E-05 
2. 144E-05 


2 

1 

Min = 2. 144299E-35 
Max = 9.02108 1E-D4 


Figure 4. Strain energy density contours near the hole boundary 

The variation of strain energy density along the boundary of the hole is plotted against angular 
position, 0 in Fig. 5. Several distinct local minima of the strain energy density function are 

observed, among which the maximum occurs at 0 = 78°. Thus, we can conclude that failure will 
initiate at the point on the hole boundary which is at a 78° angle from the horizontal axis as 
shown in Fig. 6. 

According to SED criterion the fracture propagation occurs in the direction where dW/dV reaches 
its critical value. In this case the local minima of dW/dV along various radius vectors centered at 
the point of fracture initiation need to be calculated. Thus three circles of radii ° 0.005 in, 

0.010 in, 0.015 in are drawn. Using the plot of dW/dV versus 0, as shown in Fig. 7, the maximum 
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Strain Energy Density (lb-in/in 3 ) 


of (dW/dV) min is found to occur at approximately 50°. The direction of damage propagation is 
shown in Fig. 8. 
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Figure 5. Strain energy density variation with angular position showing the location of failure 

initiation to be at 78° 



Figure 6. Point of failure initiation at the hole boundary 
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Figure 7. Strain energy density variation with angular position showing the direction of failure 

propagation to be around 50° 



Figure 8. Direction of damage propagation 

For damage tolerance consideration, a limit will be imposed on the strain energy density function to 
prevent it from exceeding the critical value associated with the material system used. 




2.2 Response Surface Methodology 

The use of finite element method for buckling analysis of a panel, representing a local model, in the 
optimization process of the global structural model would be very time consuming and inefficient. 
A more efficient alternative would be to use algebraic models derived from response surface 
methodology (RSM). 

RSM uses mathematical and statistical techniques to generate algebraic equations for a particular 
response based on the specified set of variables or factors. For example, an RS model can be 
developed for estimating the buckling strength of a panel based on a set of geometric and/or 
stiffness variables. 

Prior to applying this method to the composite fuselage structure, we applied it to the aluminum 
rectangular plate described in the previous section. The factors considered included the following 
geometric attributes: plate thickness, plate aspect ratio, and hole diameter. Since material variation 
was not considered, no stiffness parameters were included as independent factors. 

The full factorial experiment (FFE) method is used to generate the response data. The FFE 
involving k factors with each having two levels is referred to as the 2* design. With only three 
factors in this case, the FFE requires 8 (2 3 ) experiments. The experimental conditions are identified 
in Table 4. The low and high levels are denoted by -1 and 1, respectively. The low and high levels 
for these factors are as follows: plate thickness t = 0.5", 1", aspect ratio a/b =1,2, and hole diameter 
D = 0. 1 ", 0.2". The loaded edge of the plate (b) is kept constant at 10". 


Table 4. Experimental condi t ions and corresponding response values 


Experiment 

Treatment 

combinations 

Design factors 
x„ x 2 , x 3 

t, in 

a/b 

1 

D, in 

Buckling load (lb) 

1 

(1) 

-1,-1,-1 

0.5 

1 

0.1 

5.0420E5 

2 

a 

l.-l.-l 

0.5 

2 

0.1 

3.2795E6 

3 

b 

-1, 1,-1 

1 

1 

0.1 

2.4930E5 

4 

ab 

1, 1,-1 

1 

2 

0.1 

1.8129E6 

5 

c 

-1,-1, 1 

0.5 

1 

0.2 

5.0260E5 

6 

ac 

1,-1, 1 

0.5 

2 

0.2 

3.2765E6 

7 

be 

- 1 , 1, 1 

1 

1 

0.2 

2.4830E5 

8 

abc 

1,1,1 

1 

2 

0.2 

1.8108E6 


Since aspect ratio and hole diameter are geometric entities, any changes in them require the creation 
of a separate finite element model. In this case four different finite element models had to be created 
in order to account for low and high levels of aspect ratio and hole size. The finite-element based 
buckling load corresponding to each experimental condition is given in the last column of Table 4. 

The effects of x,, x 2 , and x 3 are found as 

x, = (l/4)[a + ab + ac + abc - (1) - b - c - be] = 2.1688E6 

x 2 = (l/4)[b + ab + be + abc - (1) - a - c - ac] = -0.8604E6 

x 3 = (l/4)[c + ac + be + abc - (1) - a - b - ab] = -1,925.0 

The two-factor interaction effects are found as 
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x,x 2 = ( l/4)[abc - be + ab - b - ac + c - a + ( 1 )] = -0.6058E6 
x,x 3 = (1/4)[(1) - a + b - ab - c + ac - be + abc] = -625.0 
x 2 x 3 = (1/4)[(1) + a - b - ab - c - ac + be + abc] = 375.0 
and the three-factor interaction effect is found as 
x,x 2 x 3 = (l/4)[abc -bc-ac + c-ab + b + a- (l)] = 75.0 

Examining the magnitudes of the effects shows the plate thickness (x,) to be the most dominant 
effect followed by aspect ratio (x 2 ) and the x,x 2 interaction. The effect of hole size in this case is 
rather insignificant due to its small magnitude. 

Based on the statistical analysis on the significance of each effect, a multiple linear regression 
model defined by 

P'cr = Po + P. x i + P 2 x 2 + (3i 2 x,x 2 + e (11) 

is the simplest model that could be used to estimate the plate buckling load as the response of 
interest with £ representing the random error term. In this case the random error term is zero as 
computational simulations are used in place of physical experiments to generate the response data. 
Therefore, only the first four terms are of concern in Eq. (11). The regression coefficients (3,, P 2 , 
and P| 2 are estimated by one half of the corresponding effect estimates, and (3 0 is estimated by the 
average of the buckling loads in FFE. Thus, 

P' cr = 10 6 (1.4605 + 1.0844 x, - 0.4302 x 2 - 0.3029 x,x 2 ) (12) 

is the response surface model for the plate buckling load with corresponding predicted values given 
in the fourth column of Table 5. The residual or difference between the finite-element based 
“actual” buckling loads and the estimated (fitted) values are shown in the last column of Table 5. 


Table 5. Comparison of actual and estimated axial buckling loads 


Experiment 

x 

Design factors 
x„ x 2 , x 3 

FE solution, 

Per 

FFE, P' cr 

Residual, 
e = P -P’ 

cr cr 

Standardized 
residual, d 

1 

-1,-1, -1 

5.0420E5 

5.0340E5 

800 

0.54928 

2 

1, -1,-1 

3.2795E6 

3.2780E6 

1500 

1.02990 

3 

-1, 1,-1 

2.4930E5 

2.4880E5 

500 

0.34330 

4 

l.l.-l 

1.8129E6 

1.81 185E6 

1050 

0.72093 

5 

-1,-1, 1 

5.0260E5 

5.0340E5 

-800 

-0.54928 

6 

1,-1, 1 

3.2765E6 

3.2780E5 

-1500 

-1.02990 

7 

-1,1,1 

2.4830E5 

2.4880E6 

-500 

-0.34330 

8 

1, 1, 1 

1.8108E6 

1.81185E6 

-1050 

-0.72093 


A normal probability plot of the residuals is shown in Fig. 9. The residuals appear to fall along the 
“fat pencil” line indicating the normal distribution of the residuals. 
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Residual 


Figure 9. Normal probability plot of residuals associated with the regression model 
The unbiased estimator of the residual variance a 2 is determined as 


n = 8 



n — p 8 — 4 


8.4850E6 

4 


2.12125E6 


(13) 


where SS E is the sum of squares of the residuals, n is the number of experiments, and p is the 
number of terms in the regression model. The term n - p in the denominator of Eq. (13) represents 
the degrees of freedom associated with the residuals. With the estimator of the variance known, the 
standardized residuals are computed using 



j 1) 2, w 


(14) 


with the values given in Table 5. Since all the computed standardized residuals fall in the interval of 
(-2, 2), the previously observed normal distribution of the residuals is reconfirmed. 

Although the linear regression model appears to fit the response data fairly well, it does not 
represent the physics of the problem accurately. Looking at the buckling equation for an isotropic 
rectangular plate given as 


a 


cr 


— k. 


tc^E 

12(1 -v 2 ) 



(15) 
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indicates that buckling is a nonlinear function of plate thickness, with the aspect ratio affecting the 
buckling coefficient k c . The nonlinearity of the response surface can be checked by calculating the 
response at the center point and comparing it with the average value of the responses at the low and 
high levels of thickness and aspect ratio. The finite element analysis of the plate at the center point 
corresponding to x, = x 2 = x 3 = 0 or t = 0.75", a/b = 1 .5, and D = 0. 15 gives the buckling load of 
approximately 1E6 lb. The average value for the buckling load at the four factorial points (x, = -1, 
1, x = -1, 1) is 1.46E6 lb. The difference of 46% indicates the presence of curvature at the center 
point. Therefore, in this case a second-order response surface described by the equation 

P'cr = Po H- P, X, + P 2 x 2 + P,, X, 2 + p 22 X 2 2 + p 12 x,x 2 (16) 

should be considered as an alternative to the model in Eq. (11). The central composite design 
(CCD) is used in this case to design the experiments for generating the necessary response data. In 
CCD the complete or a fraction of 2 k factorial experiments are augmented with 2k axial runs at 
value of a for each factor and at least one center point. Therefore, for 2 factors, using the full 
factorial experiments would require at least 9 experiments. The experimental conditions are shown 
in Table 6. The zero level for each factor corresponds to the average value of-1 and 1. In this case 
a is set to 2. For thickness, levels -2 and 2 correspond to 0.25" and 1.25", respectively. For aspect 
ratio, levels -2 and 2 correspond to 0.5 and 2.5, respectively. 

The estimated buckling loads at the new design points using Eq. (12) are given in the fourth column 
of Table 6. Unlike the estimates at the original design points, there are major discrepancies at most 
of the new design points. 


Table 6. Comparison of actual and estimated axial buckling loads 

Experiment 

Design factors 
x,, x„ x 3 

FE 

solution, P rr 

FFE, P'^ 

CCD, P' cr 
(1CP) 

CCD, P' cr 
(3 CP) 

CCD, P cr 
(5 CP) 

1 

-1, -1,-1 

5.0420E5 

5.0340E5 

5.022E5 

4.899E5 

4.861E5 

2 

1, -1,-1 

3.2795E6 

3.2780E6 

3.167E6 

3.155E6 

3.151E6 

3 

-1, E-l 

2.4930E5 

2.4880E5 

3.037E5 

2.915E5 

2.877E5 

4 

1, E-l 

1.8129E6 

1.81185E6 

1.757E6 

1.744E6 

1.741E6 

5 

0, 0,0 

1.0000E6 

1.46051E6 

1.059E6 

1.028E6 

1.019E6 

6 

-2, 0, 0 

4.0626E4 

-7.083 1E5 

-141.25 

5.991E3 

7.894E3 

7 

2,0,0 

4.0476E6 

3.6293E6 

4.117E6 

4.124E6 

4.125E6 

8 

0, -2, 0 

2.3146E6 

2.3209E6 

2.357E6 

2.363E6 

2.365E6 

9 

0,2,0 

7.6226E5 

6.00 14E5 

7.487E5 

7.548E5 

7.567E5 


The least squares method is used to estimate the unknown parameters in Eq. (16). The resulting 
second-order response surface equation is found to be 

P'Jicp = 1 0 6 (1-059 + 1.029 x, - 0.4022 x 2 + 0.2499 x, 2 + 0.1235 x 2 2 - 0.3029 x,x 2 ) (17) 

The estimated buckling loads using Eq. (17) are given in the fifth column of Table 6. As suggested 
by Montgomery and Runger 9 , three to five central points are usually used to develop a response 
surface model. The response surface equations using three and five central points are found to be 

P'Jjcp = 10 6 (1-028 + 1.029 x, - 0.4022 x 2 + 0.2591 x, 2 + 0.1327 x 2 2 - 0.3029 x,x 2 ) (18) 

P’ U„= 10 6 (1.019+ 1.029 x. -0.4022 x 2 + 0.2620 x, 2 + 0.1356 x 2 2 - 0.3029 x,x 2 ) (19) 

cr 5CP v 1 4 
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with the corresponding point estimates given in Table 6. Since there is no random error in 
computational experiments, the additional central points result in the same values of response as 
those shown for the corresponding fifth experiment in Table 6. It is interesting to note that the 

inclusion of additional central points affected the values of P 0 , P n , and p 22 but not the rest. 

The mean square error VISj or the unbiased estimator of the residual variance o" for each of the 
CCD-based models is determined using 


MS e = <J 2 


n 



n — p n — 6 


( 20 ) 


For one-, three-, and five-central point designs, the value of n is 9, 11, and 13, respectively. The 

corresponding a 2 is found to be 1.02285E10, 6.8039 1E9, and 4.98889E9 for one, three, and five 
central points, respectively. Using Eq. (14) the standardized residuals for each level of the CCD- 
based models are computed with the values shown in Table 7. 


Table 7. Comparison of actual and estimated axial buckling loads 


Experiment 

, ' 

Design factors 
X,, X„ X, 

FE 

solution, P,. r 

Stan, residuals 
(1CP) 

Stan, residuals 
(3 CP) 

Stan, residuals 
(5 CP) 

1 

-1, -1,-1 

5.0420E5 

0.01978 

0.17336 

0.25626 

2 

1,-1, -1 

3.2795E6 

1.11237 

1.50935 

1.81929 

3 

-1, 1,-1 

2.4930E5 

-0.53789 

-0.51160 

-0.54366 

4 

1, 1,-1 

1.8129E6 

0.55272 

0.83530 

1.01795 

5 

0, 0,0 

1.0000E6 

-0.58337 

-0.33945 

-0.26900 

6 

-2, 0, 0 

4.0626E4 

0.40309 

0.41989 

0.46342 

7 

2,0,0 

4.0476E6 

-0.68621 

-0.92622 

-1.09582 

8 

0, -2, 0 

2.3146E6 

-0.41924 

-0.58677 

-0.71356 

9 

0, 2,0 

7.6226E5 

0.13408 

0.09044 

0.07872 


The standardized residuals for all three cases fall in the interval of (-2, 2). This implies that the 
errors are normally distributed and that there are no outliers. The comparison of the standardized 
residuals corresponding to each experiment shows the quadratic response surface model based on a 
single center point to be better than the other two. However, using the root mean square error 
(MS E )° 5 as the key indicator identifies the regression with five center points (5 CP) to have the best 
fit. 


Another statistic that can be used to determine the quality of fit is the adjusted coefficient of 
multiple determination defined as 



f n - 1 "j 5S g 
K n-p) S n . 


( 21 ) 


where SS E denotes the sum of squares of the residuals, and S is the total corrected sum of squares 
of observations (P c ). Based on the data in Table 6, R 2 adj for the 1 CP, 3 CP, and 5 CP models are 
found to be 0.99484, 0.99480, and 0.99642, respectively. While there is very little difference 
between the R 2 adJ values, the 5 CP model appears to fit the data better. 
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In the coming year we plan to incorporate both the strain energy density failure criterion and the 

response surface methodology in the design of composite fuselage structure. In that case the local 

buckling load would be captured via the response surface equation and the strain energy density 

criterion will help introduce a damage tolerance constraint in the global optimization problem. 
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